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Abstract 

A two dimensional flow model is introduced with deterministic behavior consisting of 
bursts which become successively larger, with longer interburst time intervals between them. 
The system is symmetric in one variable x and there are bursts on either side of x = 0, 
separated by the presence of an invariant manifold at x = 0. In the presence of arbitrarily 
small additive noise in the x direction, the successive bursts have bounded amplitudes and 
interburst intervals. This system with noise is proposed as a model for edge localized modes 
in tokamaks. Further, the bursts can switch from positive to negative x and vice-versa. The 
probability distribution of burst heights and interburst periods is studied, as is the dependence 
of the statistics on the noise variance. The modification of this behavior as the symmetry in x 
is broken is studied, showing qualitatively similar behavior if the symmetry breaking is small 
enough. Experimental observations of a nonlinear circuit governed by the same equations are 
presented, showing good agreement. 
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1 Introduction 



This paper is motivated by observations of extreme noise sensitivity in a two-dimensional flow of 
the form 

dx 

= f(x, y) = {y - l)x, (l) 
% = g{x, y) = ty v - % 2 y- (2) 

at 

This system is a low-dimensional model for the nonlinear behavior of a plasma instability in which 
y represents the pressure gradient, and instability (with amplitude x) is driven by the pressure 
gradient and fixed magnetic field line curvature. Such pressure -driven instabilities are thought to 
be responsible for edge localized modes (ELMs) observed as fluctuations at the edge of a tokamak 
IHIZIl. Some ELMs, called Type-I ELMs, show temporal behavior which is quite simple, consisting 
of well separated large bursts, indicating that their dynamics can be represented by a low-order 
system. However, the time series appear to show chaos, and it is of some interest to determine 
whether this apparently chaotic behavior is indeed deterministic chaos or whether it is due to 
sensitivity to noise from, for example, the plasma core. For example, if the apparent chaos is 
due to noise, the behavior can occur in a two dimensional model, whereas an autonomous model 
showing similar apparently chaotic behavior must be at least three dimensional. 

The effect of noise has been studied in other experimental physics situations, and the kind 
of extreme sensitivity to noise we discuss here has been observed. For example, in experiments 
involving the formation of droplets in a viscous fluid[3|, the fluid is observed to form thin necks 
repeatedly as a part of the process. Simulations showed the formation of necks, but the repeated 
formation of necks required noise in the modeling, although extremely small noise gave agreement. 
Another example involves studies of a Nd: YAG (neodymium doped yttrium aluminum garnet) laser 
with an intercavity KTP (potassium titanyl phosphate) crystal. Theoretical studies were performed 
to model the laser dynamics[4J, showing that the type-II chaotic dynamical behavior of the laser 
was observed to be very sensitive to noise and was actually found to amplify the noise. Because of 




Figure 1: Orbits initiated near the fixed points at x = ±x = i-y/e, y — 1. The orbit on the right 
spirals out clockwise, the one on the left counter-clockwise. The x— and y— axes are, respectively, 
stable and unstable manifolds of the fixed point at the origin. 

the role of a very low level of noise in such disparate physical systems, we have been motivated to 
do detailed studies of (HI, © and related systems perturbed with a low level of noise. 

For the system ([T]), with zero noise, x grows if y > 1, but for large enough x the term —x 2 y, 
which represents the flattening of the pressure gradient due to the fluctuation, enters. This causes 
a decrease in y, which quenches the growth. For this flow, x = is an invariant manifold, and 
is in fact the unstable manifold of a fixed point at x — y — 0. See Fig. 1. The x— axis is also 
an invariant manifold, the stable manifold of the same fixed point. There are two unstable spirals 
with x = ±xo = ±y/e. The nonlinear deterministic behavior consists of spirals coming out of the 
fixed points with x = ±xo, coming closer to the two invariant axes on each pass, and developing 
increasingly larger bursts, one for each encircling of the spirals, more widely separated in time. 
Because of symmetry in x, identical bursts can occur on both sides of x = 0, isolated from each 
other by the invariant manifold x — 0. 

With a small amount of uncorrected Gaussian noise added to eq. (HJ), we find that the re- 
sulting nonlinear stochastic equation has the following property: the bursts saturate in amplitude, 
leading to behavior that is qualitatively similar to deterministic chaos. We call this behavior noise- 
stabilization. Further, the noise allows transitions across the y— axis, an invariant manifold for the 
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deterministic system. Statistically, the dynamics is symmetric. In particular, we focus on the frac- 
tion of the number of bursts with x < compared with those with x > 0; with statistical symmetry 
these are equal. In the physical system motivating this work, the processes we model as noise have 
a much shorter correlation time than the processes described by the deterministic equations ([TJ, ©, 
hence modeling them as noise is appropriate. Noise- stabilized systems are interesting for several 
reasons. Most importantly, although they can exhibit dynamical behavior that is reminiscent of 
deterministic chaos, it is likely that their behavior for very low noise level is distinguishable from 
deterministic autonomous low dimensional systems. Our model system was chosen to emphasize 
the noise-stabilizing effect, in the sense that it has no attractor in the zero noise limit. In physical 
applications, distinguishing noise-stabilized behavior from more familiar types of dynamics could 
be critical for understanding and predicting how the system under study will change as the noise 
driving is modified. 

There have been several related papers on nonlinear stochastic equations which are sensitive 
to a small amount of noise. Sigeti and Horsthemke (51 studied the effect of noise at a saddle-node 
bifurcation, and found noise induced oscillations at a characteristic frequency. Stone and Holmes 
||6l studied systems with an attracting homoclinic orbit or an attracting heteroclinic cycle (struc- 
turally stable because of the presence of a symmetry) in the presence of noise. They found that 
the effect of the noise is to prevent the time between bursts from increasing on each cycle. Stone 
and Armbruster [7| studied structurally stable (again because of symmetry) heteroclinic cycles in 
the presence of noise, and analyzed the jumping between invariant subspaces of the deterministic 
system. Armbruster and Stone JE! studied heteroclinic networks in the presence of noise, and the 
induced switching between cycles. References (BEHEl stressed the importance of the linear part of 
the flow near the saddles. Moehlis has investigated a system representing binary fluid convec- 
tion, and found that states with large bursts can be very sensitive to noise. References ifTOl ITU deal 
with a system (SEIR or susceptible-exposed-infected-recovered) describing epidemic outbreaks 
and show that chaos can be induced for parameters far from the region for which the deterministic 
system is chaotic. 
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The difference between our work and this previous work is the following. Our work concerns 
a system which, in the absence of noise, has successive bursts, each larger than its predecessor 
and separated by lengthening time intervals. In the presence of noise, our system exhibits a finite 
characteristic scale for the burst amplitude, a characteristic time for bursts, and random switching 
across an invariant manifold of the deterministic system. Further, our deterministic system is 
two-dimensional, and therefore cannot have deterministic chaos, but the noise introduces behavior 
which resembles deterministic chaos in several ways. In Refs. JUEHEl systems with homoclinic or 
heteroclinic cycles were studied; the noise was found to induce switching between subspaces and 
introduced a characteristic time scale for intervals between bursts, but the bursts in the deterministic 
system were limited in magnitude. The model of Ref. is four-dimensional, and therefore can, 
unlike our system, exhibit chaotic behavior even without noise, in principle. It was found that this 
specific system can have periodic bursts of infinite magnitude. These infinite bursts are periodic 
in the sense that if the origin and infinity are mapped to each other in a specific way, the solutions 
to the equations can reach the origin in finite time and can be integrated through it, leading to a 
periodic signal. These states with large periodic bursts were found to be sensitive to noise. This 
behavior is to be contrasted with the behavior we have found from eqs. ©, ©, in which (for v < 2) 
successive bursts get larger in magnitude, but no single burst goes to infinity, and noise causes the 
bursts to behave in a way that resembles deterministic chaos. The model in Refs. ifTOlfTTI exhibits 
noise-induced chaos because of bi-instability, related to the presence of two nearby unstable orbits. 

The model we introduce is similar to the models of Refs. J6H7J0D with a heteroclinic connec- 
tion, in the formal sense that in our model the axis is a heteroclinic orbit between the saddle at 
(x, y) = (0, 0) and the point at infinity. After a change of variables, the point at infinity can be 
mapped to a finite point and the origin can be left fixed. The new unstable manifold maps from 
the origin to this second fixed point. However, additive noise in our system would then map to 
non-additive noise in the compactified version. In particular the noise disappears at the second 
fixed point, which is physically unrealistic. 

In Sec. 2 we introduce the deterministic form of the model and show that with i/ = 1 it is 
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equivalent to the Lotka-Volterra predator-prey model. We discuss the surface of section map x — > 
x' = F(x), taking minima of x to maxima of x (and vice-versa), as well as the composite map 

x — * x". 

In Sec. 3 we introduce the stochastic model and present results. These results include those on 
the Lyapunov exponent hi and the distribution of maxima of \x\ and the time interval T between 
bursts, and the dependence of these quantities on the noise diffusion coefficient D. A brief discus- 
sion of the behavior near the y— axis is shown. In this limit, the behavior in x is linear and can be 
treated by the Fokker-Planck equation, discussed in more detail in Appendix A. 

In Sec. 4 we discuss the role of reflection symmetry in x and the effect of weak symmetry 
breaking. We also present results involving modifications to the system at small and large y, 
and a modified form of the equations in which the noise is replaced by a sinusoidal perturbation. 
The results with an offset show that in a sense the system with noise is structurally stable. The 
results with a sinusoidal perturbation lend credence to the validity of the Lyapunov exponent for 
the random case. 

In Sec. 5 we show results from an experiment with a nonlinear circuit, showing noise stabiliza- 
tion in a physical system. 

In Sec. 6 we summarize our work. 

2 Deterministic model 

The deterministic form of the model we study is eqs. CQ), 0. The parameters e, v are the only 
parameters that cannot be removed by rescaling x, y, and t. Starting with x = and y > 0, y 
increases in time, going to infinity in finite time if v > 1 . For y > 1 small initial values of x begin 
to grow. [The instantaneous growth rate of x in © equals y — 1.] If x grows at a rapid enough 
rate relative to y (to be quantified later), the second term in © eventually dominates the first and 
y decreases. For v = 1 the system ©, © is the Lotka-Volterra predator-prey model. The usual 
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form lfT2l of this system, in scaled variables, is 

as 

With the change of variables X = x 2 /2, Y = y, s = 2t, E = e/2, it can be put in the form of 
eqs. ©, © with v — 1. Notice that in this latter form there is a symmetry x — > — x not present in 
the usual form. For this value of v, equations © © can be written in terms of q = In x, p = In y 
in the form 

^ = e> - 1 
— = e — e \ 

Thus eqs. ©, © are an autonomous Hamiltonian system, with canonical variables (q,p): 

H(q, p) = e p — p + -e 2q — eq = y — In y + -x 2 — e In x. (3) 

Successive intersections of H = const, with y — 1 define a 1£) surface of section map x — > a/ = 
F(x). See Fig. 2. There are centers at x = ±x = ±y/e, y — 1. The mapping F is determined 
from H(q,p), i.e. 

-x 2 — elnx = -x' 2 — em a/. (4) 
2 2 



For small x we find x ~ x' exp(— x' 2 /2e), which can be approximated further by x' = y/—2e\nx. 
Thus for small x or very large x' , F(x) is logarithmic in nature. For large x or small x' we have 
the inverse x' = xexp(— x 2 /2e). 

On the other hand, for v > 1 the system is not Hamiltonian. It has fixed points at y = 1, 
x = ±xo= ±\/e and at x — y — 0. Near these fixed points, orbits evolve according to the 
Jacobian J(x, y) = Vf, i.e. 

^x(t) = J5x(t). (5) 
at 
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Figure 2: Contours of the Hamiltonian © in x, y for the Lotka-Volterra model \y = 1 in eq. ©], 
showing the fixed point at (x, y) = (y/e, 1) (labeled FP) and the surface of section map x — > F(x). 
For v < 1 the orbits spiral into the fixed point; for 1 < u < 2 the orbits spiral out for all time; for 
v > 2 the orbits spiral out, but as soon as they cross y with a small enough value of x they go off 
to infinity in one pass. See Sec. 2.2. 



For eqs. ©, ©, 



J(x,y) 



y — 1 x 

,V-X 



—2xy euy v — x 

For the two fixed points at x = ±x , y — 1 the eigenvalues satisfy A 2 — e{y — 1)A + 2e = 0, and 
are complex with positive real parts (unstable spirals) for 



< v - 1 < 



(6) 



Orbits continue to spiral out for v > 1. This is demonstrated by showing that the Hamiltonian for 
the case v — 1 in eq. © is a Lyapunov function for v ^ 1. To show this, we note 



dx dH dy OH 



dt 



dt dx dt dy 



(y - 1) (y^ - 1) . 



Thus, for u > 1, dH/dt > and the orbits spiral outward for all time, since H has a minimum at 
x = x ,y = 1. For u < I, dH/dt < and the orbits spiral in to the fixed point. 

The system has another fixed point, but with non-analytic behavior in y for noninteger u, at 
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Figure 3: Orbits (a) x(t), (b) y(t) and (c) phase plane y(x) for the deterministic equations <[TJ and 
©, with e = 0.5, v = 1.2, with an initial condition near the fixed point at x = y/e, y = 1. The 
orbit spirals out of the fixed point, continuing to expand, eventually piling up near the invariant 
manifolds x — 0, y = 0, with bursts to large values of x and y and long interburst time intervals 
spent mostly near x = y = 0. In (d) the finite time Lyapunov exponent hx(i) is shown. 

x = 0, y = 0. The axes x = 0, y = are invariant manifolds; we consider only y > 0, and for 
the noise-free case orbits with x(0) > remain in that quadrant. In the range of e and v given 
in eq. ©, orbits spiral away from the fixed points at (±xo, 1) [Fig. 1], approaching the x— and 
y— axes, as shown in Fig. 3, which has e = 0.5, v = 1.2. After an initial transient, the motion 
is bursty, with each successive oscillation coming closer to the axes, leading to a larger interburst 
interval, followed by a larger burst. 
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Figure 4: Zero contours of the larger eigenvalue p of the symmetrized Jacobian J s = (J + J T )/2, 
showing the fixed point (FP) (x, y) = (y/e, 1), the region (+) where p > and two regions (— ) 
where p < 0, one a very thin sliver near the x— axis. Also shown is a representative orbit spiraling 
out from the vicinity of the fixed point. The parameters are as in Fig. 3. 

We compute the finite time Lyapunov exponent 



where 8x(t) is evolved according to eq. © and x(t) is evolved by eqs. ©, <[T]). In deterministic 
systems with a chaotic attractor, hi(t) measures the average exponential rate of divergence, or 
stretching, over < t' < t. The largest Lyapunov exponent is the limit of hi(t) as t — > oo 
or the average, with suitable invariant measure, of hi(t) over the attractor. In this 2D system 
without time dependence and with diverging orbits, the infinite time Lyapunov exponent does not, 
strictly speaking, have significance. However, we will discuss hi in more detail in this section and 
Sec. 4.5, where the orbits are bounded and it is therefore appropriate. The exponent hi(t) is shown 
as a function of time in Fig. 3d. It is clear that hi [t] shows the bursts in x and y, and decreases 
whenever the orbit is near enough to the origin. In Fig. 4 we show the zero contours of the larger 
eigenvalue p(x,y) of the symmetrized Jacobian J s = (J + J T )/2, computed analytically. This 




(7) 
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quantity is relevant because |5x(t)| = (Sx.(t), 5x(i)) 



I 2 evolves according to 



(d/dt)(5x(t),Sx(t)) = (5x(t),2J s 5x(t)) < 2p{t), 



sothatp(t) = p(x(t),y(t)) is an upper bound for the local contribution to hi(t), namely (d/dt) ln|5x(t)| = 
|<&c(t)| -1 (d/dt)|£x(t)| < p(t). From this we find dhi(t)/dt < \p(t)-hi(t)]/t, (d/dt)(thi) < pit) 
or hi(t) < t^ 1 J ' p(s)ds. 

Further insight into the bursty nature can be obtained by finding the surface of section, shown 
in Fig. 5 and discussed above for the Hamiltonian case v — 1. For the parameters of Fig. 3, 
this map x — > x 1 = F(x) is shown in Fig. 5a. The slope F\x) at the fixed point x = s/e, 
computed numerically, equals si = —1.17. This value agrees with the value obtained from the 
complex eigenvalues A of J(xo,yo)> which satisfy A = A r ± i\ with A r = e{v — l)/2 and Aj = 
±V 7 Ze (1 + 0(e(u - l) 2 )), which equals ±1 for e = 0.5 and u < 1. This gives s x w -e^^" 1 ^/ 2 . 
For e = 1/2, z/ = 1.2, this gives si = —1.17, in agreement with the numerical results. This value 
si is less than —1, as it must be because the fixed point is unstable. Note that the values of x' for 
small x rise rapidly as x — > [x' is approximately proportional to V '— In x, as suggested by the 
v = 1 (Lotka-Volterra) results discussed after eq. ©], indicating that orbits that are near x = 
when they pass y = 1 lead to large succeeding maxima. Even more pronounced is that for x > 3 
the values of x' are vanishingly small, showing that moderately large maxima lead to succeeding 
minima that are extremely close to the y— axis. In Fig. 5b we show the composite surface of section 
x — * x", from one minimum to the next, or one maximum to the next. The slope at the fixed point 
is 1.37 ~ s 2 , as expected. For large x, x" = F 2 {x) appears to be exponential in x. 

Next, we turn to a discussion of the choice of the parameter v. Let us investigate the range of 
the parameters u, e for which the system exhibits successively larger, more widely separated bursts. 

Consider eqs. (HJl, for large y and small x, i.e. 



dx 



(8) 
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Figure 5: Surface of section (a) x — > x' = F(x) from one crossing of y — 1 (x — 0) to the next, 
showing x = as the fixed point. Parameters are as in Fig. 3. Composed surface of section (b) 
x — > x" = F 2 (x) = F(F(x)). The dashed lines are, respectively, x' — x and x" = x. 
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dy 
dt 



g(Q,y) » ey v . 



From these we conclude 



x = x c exp 



,2-v 



(9) 



(10) 



Le(2 - i/) 

where x c exp [l/e(2 — i/)] is the value of x when the orbit passes y = 1 with small x. Let us 
compare the two terms on the right in eq. ©, first for u = 1 (Lotka-Volterra). The second term 
exceeds the first if x 2 > e and, since x ~ e y / e , the nullcline dy/dt = is crossed, and y eventually 
decreases. For 1 < v < 2, the nullcline is crossed when x 2 > ey v ~ x or 



x 2 exp 



■2y 



2-u 



6(2 



> £2/ 



(11) 



which occurs eventually. So, in each burst, y reaches a maximum and begins to decrease, starting 
a new cycle, as long as x ^ 0. (The orbits with x = go to infinity in finite time for v > 1.) 

For = 2, we can use eq. © with eq. © for arbitrary x (including the term —x 2 y) to obtain, 
for large y, 



dy_ 

dx 



y 

e x. 

x 



The solution is 



y = (x e 



X 



with C > 0; the nullcline has y = x 2 /e. For e < 2, the nullcline is crossed and the cycle begins 
again. For e > 2 the nullcline is not crossed and the orbit can go off to infinity in one cycle, in 
finite time. 

For v > 2, the nullcline in eq. (fTTT) is never reached if x c is small enough. This means that if 
the value of x when the orbit crosses y — 1 is below some critical value, the orbit will go off to 
infinity before another cycle. Therefore, an orbit starting near the fixed point (x, y) = (y/e, 1) will 
encircle the fixed point a finite number of times and then go off to infinity in finite time. 
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3 Stochastic model and results 



3.1 Model 



With noise, the system based on eqs. (OQ), is a nonlinear stochastic ODE, of the form 



dx 
~dt 



f(x,y) + V2D£(t), 



(12) 



dy 
dt 



g(x,y), 



(13) 



with £(t) representing uncorrelated unit variance Gaussian noise, having = 0, = 

S(t — t'). Here, D is the Brownian diffusion coefficient. For a low noise level, £(t) affects the 
dynamics only near the y— axis, where f(x, y) is small. The motivation for including noise in the 
a;— equation but not in the y— equation is the following. Without noise, when the orbit is traveling 
along the y— axis for y < 1, x{t) can decrease to a level that is unrealistically small for modeling 
any physical application with noise. Noise prevents x from becoming so small for < y < 1, and 
therefore is expected to prevent the successive bursts from continuing to increase in magnitude, 
with increasing interburst time interval. We do not include noise in the ^—equation because noise 
could cause y to become negative when the orbit is near the x— axis. We will discuss a model 
allowing negative y in Sec. 4. 

We integrate the nonlinear stochastic ODE system (fTZl) . (fT3l) numerically, with a noise term in 
x added at each time step. Specifically, the time stepping from t to t + h is 



The implicit form of the deterministic part of the equations is solved by a simple Picard itera- 
tion. The random term is added after this iteration on the deterministic equations has converged. 
Each value £ (t) is an independent random number with zero mean Gaussian distribution and unit 



x(t + h) = x(t) + hf 



x(t)+x(t+h) y(t)+y(t+h) 
2 ' 2 



) 



(14) 
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variance, and the coefficient \JlDh is chosen to give results independent of the time step h (in a 
mean- square sense) in the limit h — > 0. 

3.2 Numerical results 

Results for the same parameters as in Fig. 3, with noise having D = 5 x 10~ 9 , are shown in 
Fig. 6, with < t < 1000. The orbits are still of a bursty nature, but the bursts and the interburst 
time intervals are limited in magnitude. The successive bursts appear to be uncorrelated and bursts 
with x negative are as common as those with x positive, after the transient near the fixed point at 
x = x = y 7 ?, y — 1. To the eye, these results appear similar to those of a chaotic deterministic 
system, e.g. the y — z projection of the Lorenz svstem lfT3ll . 

In Figure 7 we show the finite time Lyapunov exponent h\{t) for the case of Fig. 6 for < 
t < 10 4 . The orbits x(t) = (x(t),y(t)) given by eqs. (TTZl) . (fT3l are affected by the noise but 
the variational form for <5x(i) is eq. © and does not directly involve the noise. [Two orbits xi (t) 
and x 2 (t) = xx(t) + 5x(t) with slightly different initial conditions are integrated in time with the 
same realization of the noise £(£).] For these parameters hi(t) converges to 0.032 as t — > oo. For 
several other values of e, v, and D, with 1 < v < 2 and ©, similar results are obtained. This 
positive Lyapunov exponent shows exponential divergence between nearby orbits. This suggests 
what appears to be evident from Fig. 6, namely that the orbits behave chaotically. This conclusion 
is reasonable because the system with noise is 2D with time dependence, and because the orbits 
remain bounded for the time intervals studied, during which h\ {t) appears to converge to a constant 
value. We will return to this discussion in Sec. 4.5. 

To analyze the bursts in terms of amplitude and time interval between bursts, we introduce 
x n , x n+ \ and T n . (See Fig. 3.) These are, respectively, the amplitude (in x) of a burst (a local 
maximum for positive x, a local minimum for negative x), the amplitude of the following burst, 
and the time interval between them. In Fig. 8 we show scatter plots of T n vs. x n , x n+ \ vs. T n , and 
the composite x n+ \ vs. x n for the parameters of the case of Figs. 6 and 7, indicating the probability 
density functions f\(x n , T n ), f2(T n , x n+1 ) and f^(x n , x n+ i). These are the marginal distributions 
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Figure 6: Orbits (a) x(t), (b) y(t) and (c) phase plane y vs. x for the system with noise, eqs. (fT2t 
and (fT3"l) . The parameters are equal to those in Fig. 3, with D = 5 x 10~ 9 . The initial condition 
is near the spiraling fixed point, so that the transient spiral shows. Note that the maximum time 
t = 10 3 is much larger than in Fig. 3. 
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Figure 7: The finite time Lyapunov exponent up to t = 10 4 , for the case of Fig. 6, showing a 
positive limiting value, lim^oo hi(t) = 0.032. 

of the full distribution g(x n , T m x n+1 ) projected over x n+ i, x n , and T n , respectively. The first has 
very little scatter. This property is related to two aspects. One is the fact that the noise is added 
only to x(t) and has little effect except when x is small. The other is that most of the time interval 
T n is spent near the saddle at x = y = 0, after the burst but before the orbit can be influenced again 
by the noise, as it passes along the y— axis near y = 1. This lack of scatter shows a very strong 
correlation. However, this correlation is strongly nonlinear and would not be reflected in the linear 
correlation coefficient, but would require a diagnostic such as the conditional entropy lfT4l . The 
other plots show the expected symmetry in x. Specifically, there are four equivalent peaks in the 
four quadrants in Fig. 8c, showing that successive peaks are positive or negative, independent of 
the sign of the previous peak. Fig. 8b shows a long tail in T n , and sharp cutoffs for small \x n \ and 
small T n . 

In Fig. 9 are histograms, showing the marginal distributions of x n , at the maxima of \x\, and 
the interburst time T n . (See Fig. 3.) The maximum time was t = 10 6 and there were about 23000 
peaks in x n and the same number of interburst intervals T n . The histogram of x n is symmetric and 
shows peaks at \x n \ = 3.7, with tails around |sc n | = 4.5 and a sharp cutoff inside at \x n \ = 3.3. The 
latter histogram, reflecting the nonlinear correlation of T n with x n shown in Fig. 8a, has a strong 
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Figure 8: Scatter plots (a) T n vs. x n , (b) x n+ \ vs. T n and (c) x n +i vs. x n for the case of Fig. 6. 
Note that there is hardly any scatter in (a). The extent of the burst [measured as \x n \ or as the 
peak of y(t)] determines T n , because after a larger burst the orbit approaches the origin closer to 
the a;— axis, because most of the interburst time is spent near x — y — 0, and because the noise 
is effective only near the y— axis. The statistics plotted in (b) is symmetric in x n+1 and has a long 
tail in T n . The plot in (c) is symmetric in x n and x n+ i, with four essentially identical peaks near 
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Figure 9: Histograms Ni(x n ), N2(T n ) of (a) x n , maxima of \x\ and (b) time intervals T n , respec- 
tively, showing the marginal distributions for these quantities. The histogram of \x n \ in (a) has tail 
with \x n \ > 5 and a strong cutoff for \x n \ < 3.3; T n in (b) also has a tail to the right and a sharp 
cutoff to the left. For this case the mean values are (|a; n |) = 3.87 and (T n ) = 49.1, respectively. 
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cutoff inside T n = 30, a peak at T n = 38, and a tail for T ~ 60 - 80. 

Based on Sec. 2.2, we expect considerably different results for v > 2. These results show that, 
for the deterministic system, if the value of x at the throat y = 1 is small enough, the orbit will go 
off to infinity before another cycle occurs. Therefore, we expect that if the noise level D is small 
enough, the orbit may have a few bursts, but will diverge to infinity as soon as the cycle comes 
close enough to x = as it crosses y = 1. For large values of D, the orbit may behave as in Fig. 6 
for a very long time, but whenever x becomes small enough at the inner crossing of y = 1, the 
orbit will also go to infinity before another cycle. Numerical simulations bear this out. 

3.3 Fokker-Planck analysis near x = 

The peaks discussed in Figs. 8 and 9 are maxima in which occur at y = 1. These are related to 
the values of x near zero for which y — 1: for small values of D, the noise is important only near 
the y— axis, and as the orbit lifts off this manifold it essentially obeys the deterministic equations, 
and therefore the peaks in \x\ are determined to high accuracy by the crossing of y — 1 for small 
x. In this section we quantify this behavior by means of analysis involving the Fokker-Planck 
equation for behavior near the y— axis. 

As the orbit travels near the y— axis, x(t) satisfies the linear stochastic equation 

^ = -y(t)x + t(t), (15) 

where j(t) = y(i) — 1; for small x, y satisfies y = ey u , independent of x. The noise has 
the statistical characteristics described after eqs. (fT2l) . (fT3l) . Linearization in x holds for small D, 
up to the time when the term — x 2 y in eq. © becomes important. For low noise level (small D), 
the successive bursts are large in magnitude, leading to small values of x on the next pass. On 
each successive pass near y — 1, the correlation with the previous peak of \x\ is lost, according to 
the results shown in Fig. 8. This behavior is due to the fact that for g(0, y) = ey u with v > 1, x 
becomes small enough to become dominated by the noise while y < 1. 
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Figure 10: Sketch of deterministic orbits near x — 0. The minima of \x\ are at the throat y = 1. In 
this region, the equations can be linearized with respect to x and noise can have a large influence. 
The values y = y\, 1, y 2 correspond to t = t\, 0, t 2 hi the text. 

In Appendix A we have included an analysis based on the Fokker-Planck equation for orbits 
near x = 0, where eq. (fT3t is valid. Conclusions based on this Fokker-Planck analysis and direct 
simulations are the following. The mean value (\x n \) (c.f. Fig. 9a) decreases with D. The depen- 
dence of this quantity is shown as a function of D in Fig. 11a. The mean of the histogram of the 
interburst time T n as a function of D is shown in Fig. lib. The results for small D in Fig. 11a 
are qualitatively similar to the behavior of F(x) shown in Fig. 5a. This is expected because, as 
we have discussed in Appendix A, the orbits cross y = 1 with typical values of x proportional 
to a x ~ a 1 / 2 ~ D 1 / 2 /e 1 / 4 , and proceed with little subsequent effect of noise. The dependence 
of (\x n \) on D appears to be approximately logarithmic for small D, consistent with the approx- 
imately logarithmic behavior of the map F shown in Fig. 5a. It is also interesting to note that, 
although hi increases with D, the increase is logarithmic (for D < 5 x 1CT 5 ) and slow, varying by 
just over a factor of two for 5 x 1CT 12 < D < 5 x 10~ 4 . This logarithmic behavior extrapolates 
to hx = at the very low level D = 10~ 19 , giving \ / 2Dh = 2 x 10~ n [c.f. eq. (fT4b l. Near this 
value of D, h\ appears to begin to diverge from logarithmic behavior to remain positive. However, 
at these low noise values, roundoff is comparable to the applied noise. 

The analysis in Appendix A shows that for small x, near the intersection with y = 1, x has a 
Gaussian distribution, f(x) oc e~ x2 / 2a * . This yields a distribution for x', at the next crossing of 
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Figure 11: Mean values of (a) the burst peaks (|x n |), (b) the interburst time (T n ) and (c) the 
Lyapunov exponent hi as functions of D. The parameters (except for noise level) are the same as 
in Fig. 6. The quantities < \x n \ > and hi appear to behave logarithmically for small D. 
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y — 1 where \x n \ is a maximum, equal to 



g(x') = \dx/dx'\f(x(x')), 



where the functional form for x(x') is shown in Fig. 5a. The second factor is responsible for 
the sharp cutoff to the left of the peak in x' (Fig. 9a), corresponding to x being in the tail of the 
Gaussian. The tail to the right of the peak in Fig. 9a is due to the Jacobian factor \dx/dx'\. For 
example, for v — 1.2 the behavior for small x from Fig. 5 is similar to that for v — 1, derived after 
eq. ©, namely x' ~ >/— In x. From the Gaussian form for f(x) we obtain \dx/dx'\ ~ x'e _a: ' 2 and 



The first (Jacobian) factor x'e - ^ gives a Gaussian-like tail for large x' and the second factor gives 
a cutoff for x' close to the fixed point x' = x = y/e, where x' — x = —si (x — x ). This cutoff is 
sharp if cr x < i . 

4 The role of symmetry and relation with other models 

We have commented that the system (fT2"l) . <fi"3l) has certain features that are not generic. These 
issues are (a) the reflection symmetry of the equations in x; (b) the fact that deterministic orbits 
eventually go to infinity, and (c) the non-analytic behavior of y v near y = 0. In this section 
we discuss results obtained when the system is modified in these areas. To deal with issue (a), 
we destroy the symmetry in x by an offset [a constant term added to eq. (fT2l> l. These results 
suggest a modification to the notion of structural stability in the presence of noise: the behavior 
is qualitatively unchanged if the offset is small relative to the noise. To deal with issue (b), we 
show results in which the behavior for large y is modified, preventing orbits from going to large y. 
Regarding issue (c), we modify the system near y = to remove the non-analytic behavior there. 
We also discuss modifications breaking the reflection symmetry in x together with limiting the 
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behavior for large y. Finally, we discuss modifications to the system involving adding a sinusoidal 
perturbation to eq. (HJ) in place of the noise term. In these studies, conventional deterministic chaos, 
characterized by a Lyapunov exponent, is observed and compared with results with noise. 

4.1 Breaking of the symmetry in x 

We have investigated the effect of breaking the reflection symmetry x — > — x in eqs. (fT2l . (fT3"l) . mo- 
tivated by the experimental results shown in Sec. 5.3. The simplest way of breaking this symmetry 
is to introduce a constant offset. With this offset, eq. (fT2l takes the form 



with the y— equation unchanged. Numerical results with zero noise show that for a > a stable 
limit cycle is formed to the right of x = 0, and points near (x, y) = (0, 0) go into this limit cycle. 
(For a < the results are identical, with x — > —x.) Therefore the zero noise results of Sec. 2 are 
not structurally stable with respect to such an offset. 

However, in the presence of noise, the results change considerably. In Figs. 12a,b we show 
x(t) and the phase portrait y vs x for a case with the same parameters as in Fig. 6 (in particular 
with D = 5 x 10~ 9 ), but with a = 5 x 10~ 5 . The results are qualitatively similar to those in 
Fig. 6 except that most of the bursts go to the right. In Fig. 12c we show the fraction $ of bursts 
that go to the left as a function of the offset a for three values of D, and in Fig. 12d we show the 
Lyapunov exponent h\. For a < VD, hi and the fraction $ are appreciable and the orbits behave 
qualitatively as in Fig. 6. For a > \/D, on the other hand, virtually all the orbits go to the right 
($ ps 0) and have negative Lyapunov exponent and therefore behave qualitatively as the limit cycle 
found for D = 0, a > 0. These results, and those of Appendix A showing o x ~ VD, indicate that 
the offset changes the results qualitatively if it moves the orbit outside the region near x = where 
noise dominates. 

This brings up the issue of structural stability of the behavior observed for a = 0. For zero 
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Figure 12: Results with an offset [c.f. a in eq. (fT6l) l. In (a), (b) are x(t) and y vs x for parameters 
as in Fig. 6, again with D = 5 x 1CT 9 . In (c), (d) are the fraction $ of bursts to the left and the 
Lyapunov exponent h±, for three values D = 5 x 10~ 9 , 5 x 10~ 7 , 5 x 10~ 5 . 
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noise, this behavior, seen in Fig. 3, is certainly not structurally stable. However, for D > the 
qualitative behavior persists as long as a < \f~D. In this modified sense, the system with finite 
noise is structurally stable. 

We will return to the issue of an offset in the electronic circuit in the next section. 



4.2 Modifications for large y 

We have discussed the deterministic model for v > 2 in Sec. 2, showing that orbits go to infinity 
after a few passes near the fixed point (x, y) = (xq = y/e, 1). The dynamics in the presence of 
noise is the following: if the noise is large enough, the value of x at the throat where y = 1 will 
typically be large enough that the system encircles (xo, 1) many times. Even with noise, however, 
eventually an orbit comes through the throat with small enough x for the system to go to infinity 
before another cycle can occur. 

A system related to eqs. ([T]), © with orbits that do not to to infinity is the predator-prey system 
of Odell lfT5lfT2ll . This system can be put into the form 

dY 

?± = Y 2 (l-Y)-XY, 
as 

or by a change of variables (X = r]x 2 /2, Y = -qy, s — 2t/rf) 

5 = <*-!)*, 07) 



f t = «/"(! " TO) - 1 2 », (18) 
with e = v — 2, i.e. the form of eqs. (QJ, with u = 2 and y 2 — > y 2 {\ — rjy). This system 



has fixed points at x = ±y e(l — rj), y = 1. For e = v = 2 these fixed points are unstable if 
rj < 1/2 and oscillating (complex eigenvalues) if r] < y/3/2. This system also has a saddle at 
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x = y — 0, with zero eigenvalue in the y direction. In addition, it has a fourth fixed point, with 
x = and y = I/77. This fixed point is a saddle, stable in the y— direction and unstable in the 
x— direction; the section of the y— axis with < y < I/77 is a heteroclinic line. Because of the 
presence of this saddle, there are two stable limit cycles, related by the reflection symmetry in x, 
to which typical orbits converge. For r] small, this limit cycle has large excursions, with peaks in 
y approaching 1/77. We have studied eqs. (fT71) . (fTSb with noise in x, and with e, v in the range of 
parameters of Fig. 3. The results are similar to those of (fT2l . (fT3l) . as long as D is large enough 
that the excursions almost always have y <C I/77. Specifically, the value of hi and the probability 
density plots as in Figs. 8, 9 are essentially identical. The effect of positive 77 is similar to the effect 
of clipping the voltage corresponding to y in the circuit (see Appendix B), except that by design 
the clipping turns on much more rapidly than the factor (1 — rjy) in eq. (fT8l . 

4.3 Modifications near y = 

We have studied the system (fT2l . (fl"3t with ey u — > go(y) = e(f3y + y u ). This modification regular- 
izes the vicinity of y = 0: the saddle at the origin is no longer dominated by y u , and has eigenvalues 
— 1, ef3. The spiraling fixed points have x = ±20 = ±-y/e(l + /3), y — 1. We have found that 
noise has the same qualitative influence for positive (3 as it does for (3 = 0. In Fig. 13a we show the 
scatter plot x n x n+ \ for D = 5 x 1CT 5 , e = 1.5, v = 1.2, and /3 = 0, (3 = 1 superimposed. For 
(3 = 1 the eigenvalue e(3 > 1, which implies that, when following a deterministic orbit along the 
x— axis and up along the y— axis, it ends up further from the y— axis than it started from the x— 
axis. (For the equations linearized about the origin, x ef3 y is constant.) This is related to the liftoff 
phenomenon of Refs.0|H|. Based on this consideration, one might expect that the sign of x n+ ± 
might correlate with the sign of x n , and the symmetry of the scatter plot would be broken, with the 
distribution fs(x n , x n+ i) having more points in the NE and SW quadrants and fewer in the SE and 
NW quadrants, while of course still preserving the symmetry in the marginal distribution of x n , 
J fs(x n , x n+ i)dx n+ i. Nevertheless, the scatter plot for (3 = 1 appears to have the same symmetry 
as for (3 = 0. 
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Figure 13: Scatter plot (a) x n — > x n +i in me notation of Fig. 3 (successive maxima of for 
e = 1.5, v = 1.2, D = 5 x 1CT 5 and both j3 = and (3 — 1, showing four-fold symmetry 
in both cases. Surface of section x — > x" = F 2 (x n ) (b) for < x < ye(l + (3) = 1.22 
for the deterministic case with e = 1.5 and both values of j3. For (3 = 1 the fixed point is at 

x = y/e(l + P) = 1.73. 
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This four- fold symmetry is explained by Fig. 13b, which shows the surface of section x — > 
x" = F 2 (x) for < x < xq, similar to that in Fig. 5b. For both cases x" <C x. For these 
parameters x" ~ x 3 for j3 = 1, while x" goes to zero faster than any power when j3 = 0. The 
origin is so very attracting for F 2 because small x maps to large x' under F and the orbit from there 
passes extremely close to y — 0, thereby leading to extremely small x" in spite of e(3 > 1. Because 
of this property, if an orbit starts with x ~ o x at y = 1 and executes one cycle, the value x = x" 
when it crosses y = 1 after this cycle will be so small (x" ~ cr 3 for /3 = 1) that it is dominated 
by the noise added for small x and will even for (3 — 1 be nearly independent of x. This four-fold 
symmetry was observed for these parameters for 5 x 10~ 13 < D < 5 x 1CT 3 . 

Although the increase of (3 has no effect on the symmetry of the scatter plot x n — > x n+ i, it has 
a profound influence on the burst intervals T n . For larger j3, typical values of T n (not shown) are 
much smaller because of the liftoff phenomenon. This dependence of T on (3 is understood easily. 
Suppose the orbit enters the region [0, a] x [0, b] with y = y . We find that if f3y y v , the time 
to exit the region equals T\ = (1 / e/3) \n(b/y ) ~ f3~ l . If, on the other hand, v > 1 and the orbit 
is far enough from the origin that y u ^> f3y, then (3 can be neglected and the time interval equals 
T 2 = (yo 1 ^ 1 ' - fr 1 " 1 ') / [e(i/ - 1)]. For example, for e = 1.5, (3 = 1, i/ = 1.2, 6 = 1, j/ = 10" 4 we 
find! 1 ! = 6.1 andT 2 = 18. 

We have considered other models in which go(y) is linear in y near y = but behaves as ey u for 
large y. The cases investigated were go{y) = cy(P p + y p( - v ~ l >) l /' p for various values of p, including 
p — 1. [Note that (yf is analytic at y = if p(z/ — 1) is an integer.] The results for all the tested 
values of p are similar to the p — 1 case described above. 

We have also considered the case go(y) = e((3y + y 2 ). In Sec. 2 we concluded that the de- 
terministic system for v — 2 continued to have bursts of increasing amplitude and time interval 
(rather than being capable of going to infinity in finite time in a single burst) if e < 2. Results 
for various values of e < 2 and D show that the results are similar to those for v < 2, as long 
as e is small enough, [3 is large enough, and D is large enough. Note that for this case the flow 
is analytic everywhere, including y < 0, and that for the deterministic form there is a fixed point 
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at x = 0, y = —(3, as well as the fixed point at the origin, the latter having the a;— axis as its 
stable manifold. This new fixed point is attracting in both directions, and therefore any noise in the 
y— direction eventually leads the orbit to this fixed point. 



4.4 Limitation for large y with asymmetry in x 

A phase portrait for a flow including both saturation in y and symmetry breaking in x [c.f. eqs. (fTBTl . 



is shown in Figure 14, with 77 = 0.1 and a = 0.015. The unstable spirals are now slightly asym- 
metric due to the finite value of a. There are two saddle points near the origin, at x = a, y = and 
at x ~ a, y « (a 2 / e) 1 ^ v ~ 1 \ For these parameters, the second fixed point has y ~ a 10 /e 5 ~ 10~ 17 , 
and the dynamics of the system can be described as if only the saddle at x — a, y — exists. 
This saddle still has the x— axis as its stable manifold (with right and left pieces labeled 1SR and 
1SL). The unstable manifold of this saddle (labeled 1U) now is no longer the y— axis, but bends 
slightly to the right and eventually asymptotes to a limit cycle (not shown) orbiting the right spiral 
fixed point. Another saddle at approximately x ~ — ar), y ~ I/77 (filled circle) has an unstable 
manifold with right and left pieces (labeled 2U R and 2UL, respectively). The invariant manifolds 
bend downward, coming into the vicinity of the x— axis, pass very close to the saddle at the origin, 
and both converge onto the unstable manifold 1U, thus approaching the limit cycle on the right as 
well. The stable manifold for the upper saddle point, labeled 2S, if followed backward in time, 
asymptotes to the spiral on the left. Hence, a narrow region on the y— axis near y = 1 that is 
bounded by 2S on the left and 1U on the right sets the scale for the noise response. If the noise 
amplitude is smaller than the width of this region (denoted A), nearly all points passing through 
this region will go to the right and asymptote to the limit cycle. If a x > A, then orbits will get 
kicked to the left and right with nearly equal probability, leading to noise stabilized behavior that 
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Figure 14: Phase portrait for the deterministic system of equations with limiting in y and offset in 
x, with zero noise and v = 1.2, e = 0.5, a = 0.015, 77 = 0.1. There is a saddle near the origin 
(open circle), a saddle near y = 1/rj (filled circle), two unstable spirals, and a stable limit cycle 
(not shown) on the right. The symbols 1SL, 1SR, and 1U represent the left and right arms of the 
stable manifold of the fixed point near the origin and its unstable manifold. The lower arm of the 
stable manifold of the fixed point near y = 1/r/ is 2S and its unstable manifold is 2UL, 2UR. 
Points from the spiral on the right go to the stable limit cycle; points from the fixed point on the 
left eventually end up outside 2S and go to the same limit cycle. 
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prevents the relaxation onto the limit cycle. 

Thus, the presence of the symmetry breaking term in the deterministic dx/dt equation destroys 
the heteroclinic connection between the two saddle points, leading generically to a limit cycle 
either on the right or left, depending upon the sign of the offset. Thus, the deterministic dynamics 
for a = 0, 77 > discussed in Sec. 4.2 is not structurally stable, but the behavior with noise is 
structurally stable in the sense discussed at the end of Sec. 4.1. The noise response is very similar 
to the noise response of the model with rj = a = (Sec. 3), to the model with 77 = 0, a ^ (Sec. 
4.1) and to the model with 77 > 0, a = (Sec. 4.2). 

4.5 Sinusoidal perturbation 

We have integrated eqs.©, © with a sinusoidal term £(£) = bsm(ujt) added to the x— equation 
rather than random noise. We chose uj to be large enough so that the sine goes through many cycles 
when the orbit is along the a;— axis, but large enough to avoid aliasing, i.e. uh < ir, where h is 
the time step. The sinusoidal and random forms of are extremes of temporal driving, with 
quasiperiodic time dependence and colored random time dependence as intermediate cases. In all 
such cases the analysis of Sec. 3.2 indicates that the typical value of x at y = y 2 is the important 
factor. (See Sec. 3.2 and Fig. 10.) This suggests that the Lyapunov exponent h\ has validity in all 
these cases. To explore this further, we have obtained results for v = 1.2, e = 0.5, as in Fig. 6, 
and with various values of u and b. The results were found to be qualitatively similar to those 
with noise, with a simple relation between b and D, showing that indeed the accumulated effect 
on x at the time y = y 2 is the determining factor. That is, o x ~ b/uj or b/uo ~ D 1 ^ 2 /e 1 / 4 . In 
particular, the behavior of < \x n \ >, < T n > and hi are similar. Thus, the similarity of the results 
with this deterministic non-autonomous system and the nonlinear stochastic system (fT2l . (fT3t lend 
credence to the idea that hi as defined in Sec. 2 and used in Sec. 3.1 is the appropriate form of 
the Lyapunov exponent for the stochastic system. It is known that a system with periodic driv- 
ing can be distinguished from an autonomous system or one with more complex temporal driving 
by means of nonlinear symbolic time series analysis! 16|. This distinction is possible because of 
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definite dips in the conditional entropy of symbolic time series when the sampling time equals 
the period 27t /l<j ||T61 . This condition distinguishes periodic driving from all other temporal driv- 
ing (autonomous, quasi-periodic, colored noise, white noise), but does not distinguish the other 
possible varieties from each other. This topic is outside the scope of the present investigation. 

5 Electronic circuit 

In order to test for noise stabilization in a physical system, we have constructed a circuit which 
integrates eqs. (13) and (14). In dimensionless integral form, these equations are x(r) = xq + 
It ((y ~ l) x £(, T ')^dT' and y(r) — yo + £^ (ey u — x 2 y)dr' ', and the parameter values used in 
the circuit were e = 0.5 and v = 1.2, as in Figs. 3,5-9,1 1,12. The circuit design is shown in Fig. 18. 
The white noise, £(£) = \/2D£(t), stabilized the oscillations, and Figs. 15-17 show that the circuit 
output agreed well with numerical solution of eqs. (13) and (14). We also observed the structural 
instability in these equations. See Appendix B for a description of the circuit design. 

5.1 Properties of the added noise 

The noise was generated by creating random numbers and recording them to a . wav file to play 
back via the computer's audio output at the standard rate of 44 kHz. This net process effectively 
filters the noise through a lowpass filter. When we sampled the noise using a digital oscilloscope, 
we found that the noise had a relatively constant spectrum to frequencies as high as 20 kHz. We 
autocorrelated the noise, and found that it was well represented by: 



with a period T = 50 /is, which also represents a flat spectrum filtered by a 20 kHz low-pass filter. 
For times longer than T/(2n), this autocorrelation function is a good approximation of A S(t). 
By evaluating the autocorrelation function at t = 0, we can determine that A = j (V^) so the 
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Figure 15: Circuit output (dots) compared to numerical solution of the ODE (traces), with pa- 
rameters as in Fig. 3. Adjusting the simulation parameters to fit the data showed that all circuit 
parameters are within 3% of their expected values. The insets show the agreement of the (digitized) 
data and simulation near the fixed point. 
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in terms of the scaled variables used in Appendix B. The theoretical minimum diffusion constant 
for our circuit parameters given by eq. (IB- 111 is well below the intrinsic noise in the circuit. This 
intrinsic noise is not well characterized and occurs in both the x and y variables. We use a large 
enough value of the noise amplitude so that the intrinsic noise contribution is negligible. We show 
in Figs. 16 and 17 the quantities T n _i vs x n and T n vs x n , first obtained from the experiment and 
also by integrating numerically the differential equations with the same parameters, in particular 
D = 4.7 x 1CT 4 . (These results are similar to those in Fig. 8, but with a different value of D.) The 
agreement is very good. 
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Figure 16: Comparison of peak height x n to (a) time since previous peak T n ^\ and (b) time until 
next peak T n , from experiments. The correlations seen here are indicative of noise stabilization. 
The noise level is D ~ 4.7 x 10~ 4 . 
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Figure 17: The same quantities as in Fig. 16 from numerical computation of eqs. (fT2l) . (fT3l . 
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5.2 Offsets and symmetry breaking 

The primary difficulty in designing this circuit is that small DC offsets at the input of the integrators 
significantly change the differential equations. In particular, an offset in the input to the y-integrator 
either drives the V^-output negative to create an error in the AD538 computational unit, or it leads to 
a stable limit cycle similar to that described in Sec. 4.4. We adjusted a small current (~ 0.45/xA) to 
minimize the V y -offset, using the automatic reset circuit to recover whenever V y became negative. 
The reset kicks the circuit back into the vicinity of one of the unstable spirals. The x-integrator 
naturally follows, bringing V x to a value near its fixed point. Without this reset, a negative value of 
V y leading to the failure of the AD538 causes the circuit to fall to a stable fixed point with a large 
negative value of V y . An external trigger can also reset the circuit to values near its unstable fixed 
point. 

Similarly, we also corrected the offset in the x-integrator by adding ~ 0.2 /xA at the integrator 
input. We adjusted this value until the noise signal generated equal numbers of negative and 
positive x pulses. After these adjustments, we observed the basic structure of the oscillations as 
they evolved away from the fixed point, in order to verify that the circuit waveforms were the same 
as the model calculations (see Fig. 15). The fact that such a simple adjustment can give results in 
agreement with the symmetric model is consistent with the extended concept of structural stability 
discussed at the end of Sec. 4.1. The results also show that the circuit is a sensitive detector of 
offsets. 

6 Summary 

We have performed a study of a nonlinear stochastic ODE whose deterministic form has unstable 
spirals, leading to bursty behavior, with successive bursts growing in magnitude and with larger 
time intervals between them. This bursty behavior is due to the fact that after each burst, the 
orbit comes closer to the unstable manifold (y— axis) of a hyperbolic fixed point at the origin, and 
therefore travels farther along this unstable manifold before diverging from it to form the next 
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burst. 

In the presence of noise at a very small level, the bursts get stabilized in the sense of becoming 
limited in magnitude. The time interval between them also limited, and the bursts can go to either 
positive or negative x. In many qualitative senses, the behavior appears like deterministic chaos. 

This system has reflection symmetry in x; an offset a in a; destroying this symmetry can lead to 
completely different behavior, depending on its magnitude relative to the noise. That is, the bursty 
behavior seen in the symmetric deterministic equations is not structurally stable. With noise and a 
small value of the offset \a\ < \JlD (D is the Brownian diffusion coefficient), the bounded bursty 
behavior persists, but with more bursts going to the right if a > (to the left if a < 0.) For larger 
offset a > \/2D, all bursts go to the right and basically give a noisy form of the stable limit cycle. 
In this sense, the results in the presence of noise and a = are structurally stable. 

We have considered modifications to the model allowing for saturation of y, because bursts 
cannot continue to grow without bound in a physical system. We have also considered modifica- 
tions near the saddle at the origin, to give the saddle at the origin a positive eigenvalue. This change 
in the linear part of the flow near the saddle affects the time intervals between bursts, making their 
characteristic value much smaller, but does not affect the properties of the burst amplitudes, or the 
signs (in x) of the bursts. 

We have described briefly results on a nonlinear circuit satisfying the same equations as the 
model. The circuit behaves similarly to the model. In particular, the circuit is very sensitive to 
the presence of an offset, and in practice the offset is adjusted to minimize the asymmetry of the 
signal. More details are presented in Ref. ifTTI and in Appendix B. 

The system (TT21) . (fT3l and its generalizations in Sec. 4 are arguably the simplest realizations 
of systems in which a small noise level can limit the amplitude of bursts and lead to qualitatively 
distinct behavior. We have listed in the Introduction physical examples of systems in which this 
effect may be important. For the tokamak example, the results here should have an impact on low 
dimensional modeling of ELMs. Indeed, the observation of chaotic time dependence of ELM data 
suggests that a simple autonomous ODE model must be three-dimensional. However, tokamaks 
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are known to have a broad spectrum of fluctuations (turbulence). If these fluctuations can be 
treated as uncorrected noise, i.e. if their correlation time is much shorter that ELM time scales, it 
is justifiable to explore two-dimensional models with noise such as the models studied here. 



Appendix A: Fokker-Planck Equation 



The stochastic behavior of eq. < fl3T > is governed by the Fokker-Planck equation for the proba- 
bility density function f(x, t), 

d f ■ 6 Mt)xf)= °( D ?f), (A-l) 



dt dx dx \ dx r 

where D = a 2 /2 is the diffusion coefficient. For arbitrary j(t), (IA-11) has the exact solution 



f( x t) = J - e -^ 2 /2«W 



if the variance or temperature a(t) satisfies 

d 



d a{t)=2[a{t) 1 {t) + D}. (A-2) 



Eq. (IA-21) has the solution 

a[t) = 2D f d Sl e 2 ^ l{s2)ds \ 

assuming a(t — > — oo) = 0. Thus, a(t) is proportional to D, with a coefficient depending on "f(t). 
If 7 is approximately constant (|7/7 2 | 1) and negative, a approaches a slowly varying state 



with a(t) = D/Y){t)\, in which the inward motion due to the advective term in (IA-11) balances 
diffusion and df /dt is negligible. This limit gives 

f{x, t) - y/WfaDe-™*' 213 . (A-3) 
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Another limit is recovered by neglecting ^(t) in eq. (IA-11) . giving 



a(t) = a(ti) + 2D(t - = 2Dt + 2a , 
where t\ is the time when 7 ~ 7 2 . Without loss of generality we can set the time where 7 = 



to t = 0. This range, in which the advective term in eq. (IA-11) is small, gives the purely diffusive 
random walk result 

f(x, t) ~ 1 -xV4(^oo). ( A-4) 

y/4n{Dt + a ) 

A third range has 7 positive with advection dominating diffusion. We find 

a(t) = a(t 2 )exp ^2 jf 7(s)dsJ , (A-5) 

where £2 is the time this range is entered, i.e. where 7(^2)^(^2) ~ D. In this range the noise 
becomes negligible. 

A simple example having these properties has 7 linear in time, j(t) = 'jot. Again taking 
a(t = —00) = 0, we find 

a(t) = 2De^ j e'^ds. 
J —00 

In this example a(t) has slow growth for t < t± = — l/V7o, diffusive increase for t± < t < t 2 , 
where t 2 — 1/ V7o, and exponential growth for t > t 2 . The value of a(t) at t — (corresponding 
to y = 1) is a 2 = a(0) ~ .0/^- 

For application to eqs. (fT2l . <TT~3b . consider x small so that its equation is linear (when the 
second term on the right in (fT51) is negligible). We then note that if a is small for y « 0, then 
a(£) near y = 1 (recall 7(t) = y(t) — 1 = 0) is proportional to D/y/jo. Since 7 = y ~ e, we 
have a(y w 1) ~ D/y/e. After a diffusive stage, a continues to increase as in eq. (IA-51) . with 
noise no longer playing a role. Thus, the nonlinear orbit for later times depends only on the noise 
accumulated by the time (here t = t 2 ) just after the orbits cross the throat at y — 1; the value of x 
at V ~ IJ2, when noise last plays a role, is proportional to \fa oc D 1 / 2 /e 1 / 4 . See Fig. 10. Thus, in 
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essence, the orbit from the crossing of y — 1 with small x out to the next crossing and back to near 
the origin is deterministic, and the noise plays its role only along the y— axis. 



Appendix B: Circuit Design 



The design of our circuit is basically the same as reported in Ref. |[T71 . but we have adjusted our 
circuit parameters, and extended the analysis of the circuit behavior. For the sake of completeness, 
we have included all of the new circuit parameters in this appendix, as well as our analysis of the 
minimum noise amplitude necessary to keep the circuit from saturating the circuit elements. 

The analog circuit consists of three basic sub-circuits: the ^-integrator, the y-integrator, and 
the reset controller, as shown in Fig. 18. The integrators use OPA 4228 operational amplifiers (low 
noise, 33 MHz bandwidth) with capacitive feedback (10 nF) to integrate their inputs. V\ and V 2 
are constant applied voltages, while V x and V y are time varying voltages, proportional to x(r) and 
y{r), respectively. 

The input to the y-integrator uses an AD538 real-time computational unit (400 kHz bandwidth) 
to raise the V y voltage to a fractional power, V^t)" -1 , by taking its logarithm, scaling the result 
by v — 1, and then exponentiating to generate Vi(V y {t) /V 2 ) v ~ 1 . This output is then added into the 
output of an MPY634 precision multiplier (10 MHz bandwidth) that creates the ratio V£(t) jVi. A 
second MPY634 multiplies this combined signal by VyjVi before it enters the integrator. We also 
use additional small adjustable current sources to eliminate offsets. 

The input to the ^-integrator is the sum of V x , the noise source, and V x V y /V2 formed by another 
MPY634. The net output signal of the entire circuit has a maximum frequency of 2 KHz, well 
within the bandwidth limit of all the components. This circuit does the following integrations: 

14(f) = VM + Si («gf> - 1) vM-jfa + SI iMOslv 
v,(t) = V4W + II (v, (qp)" - V„( ( <)) JL., 
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Figure 18: Circuit diagram. 
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where the circuit components had the values listed in Table 1, and the parameter v — 1 was set 
to 0.2 in the AD538 component by a voltage divider composed of a 2200 f2 resistor and a 560 f2 
resistor. This dimensional form of the equations is related to the dimensionless form by defining 

x, y, t, e and r\ as: 



,. _ Ra Yz 

y r 3 v 2 



T 



R2C2 

R2C2V1 I R3-' 

RiCi V 2 \ R 2 



R2C2 V x _ r- V x / R.2 



~ — / R2C2 V x _ r- 



RiCx V 2 v VWV2 \R3 



2 



R2C2 Vn R2 
RiCi V 2 Ra 



This leads to fixed points at: 



Vx* = v / VV^(f 



v-l 

2 



Thus, a circuit design with a given value of e has its fixed points and its voltage scaling determined 
by the choice of the ratio R3/R2. This value can be optimally set by forcing both the x circuit 
and the y circuit to reach saturation values on the same cycle. For the v = 1 case, neglecting the 
logarithmic terms of the Hamiltonian H(x, y) in eq. ©, the peak value of y (y p ) and its following 
peak value of x (x p ) are related by x 2 p = 2y p if H is large enough, i.e. for bursts with x p , y p large 
enough. These two peak values cannot require voltages in excess of V 2 , or the multipliers will fail, 
and the peaks will be clipped. To optimize, we equate these peaks when they reach V2; for the 
v — 1 case this gives eV^/Vi = 2R2/R3 or, for our values of V\ and V2, 

R 2 e V 2 



This choice then implies maximum values of x m = y2 (eV 2 /2Vi) = 3.53, and y m = eV 2 /2Vi = 
6.25. These maximum values of x and y determine the minimum noise amplitude that must be 



43 



Vi OAV 

V 2 10V 

R x 6.8ktt 

R 2 122fcO 
R 3 

R 4 67ktt 

C x lOnF 

C 2 lOnF 

Table 1: Values of circuit elements. 

present to keep the voltage peaks within the operating range of the multipliers. The logarithmic 
dependence observed in Fig. 1 1 can be approximated as (x) = (1/8) In (10 5 /D), so that: 

D min = I0 5 e~ 8x ™ = 10 5 e"V 2 (^) ~ 2 x 10" 10 . (B-l) 



When the amplitudes are low enough to avoid clipping, the measured results are in agreement with 
those given in Sec. 3.2. 
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